function F_sum = convolution_cdf(z, cdf_F1, c1, c2, Gamma)
    % convolution_cdf - Calculates CDF of X + Y - Gamma (vectorized using arrayfun)
    %
    % Inputs:
    %   z - point(s) at which to evaluate the CDF (can be a vector)
    %   cdf_F1 - function handle to the user-defined CDF
    %   c1, c2 - lower and upper bounds of the distribution
    %   Gamma - constant to subtract (must be < 2*c1)
    %
    % Output:
    %   F_sum - CDF value(s) at point(s) z

   

    % Ensure gamma < 2*c1
    if Gamma >= 2*c1
        error('Gamma must be less than 2*c1');
    end
    
    % Define the integration function
    f_c            = makedist("Exponential", "mu", 40);  % Cost distribution (Exponential with mean 40)
    f_c            = truncate(f_c,10,100);               % truncate c>10. 
    f_pdf          = @(c) pdf(f_c, c);

    function F = integrate_cdf(zi)
        integrand = @(y) cdf_F1(zi + Gamma - y) .* f_pdf(y);
        F = integral(integrand, c1, c2);
    end

    % Use arrayfun to apply the integration to each element of z
    F_sum = arrayfun(@integrate_cdf, z);
end